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ABSTRACT 

We show that the motions of supergranules are consistent with a model in which they are 
simply advected by the axisymmetric flows in the Sun's surface shear layer. We produce a 10- 
day series of simulated Doppler images at a 15-minute cadence that reproduces most spatial 
and temporal characteristics seen in the SOHO/MDI Doppler data. Our simulated data have 
a spectrum of cellular flows with just two components - a granule component that peaks at 
spherical wavenumbers of about 4000 and a supergranule component that peaks at wavenumbers 
of about 110. We include the advection of these cellular components by the axisymmetric flows - 
differential rotation and meridional flow - whose variations with latitude and depth (wavenumber) 
are consistent with observations. We mimic the evolution of the cellular pattern by introducing 
random variations to the phases of the spectral components at rates that reproduce the levels of 
cross-correlation as functions of time and latitude. Our simulated data do not include any wave- 
like characteristics for the supergranules yet can reproduce the rotation characteristics previously 
attributed to wave-like behavior. We find rotation rates which appear faster than the actual 
rotation rates and attribute this to the projection effects. We find that the measured meridional 
flow does accurately represent the actual flow and that the observations indicate poleward flow 
to 65° — 70° latitude with equatorward counter cells in the polar regions. 

Subject headings: Sun: granulation, Sun: rotation, Sun: photosphere, Sun: surface magnetism 

Introduction one day, and flow veloc ities of 300 m s _1 

(|Rieutord fc Rinconl l2010h . They cover the en- 



Supergranules are cellular flow structures ob- tire surface of the gun except for the immediate 

served in the solar photosphere with typical di- surroundings of sunspots 



ameters of about 30 Mm, lifetimes of about „ , ,. , , — ~n „ r .k 

Supergranules were discovered by IHart I ([1954) 
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but it was iLeighton et al. I (|1962r ) who suggested 
that these cellular flows were convective structures 
and coin e d the term "supergranule." However, 
Worden (1975). and others since, had difficulty 
detecting any associated thermal features (i.e. hot 
cell centers) consistent wit h a thermal con v ection 
origin. Recently however, Meunier et al. I ( 2007 ) 
have reported the detection of a small temperature 
excess at cell centers of about 0.8-2.8 K. 

Supergranules are intimately involved with the 
structure and ev olution of the magnetic fi eld in 
the photosphere. ISimon fc Leighton I (|l964l ) found 
that the magnetic structures of the chromospheric 
netw ork are located at the boundaries of these 
cells. Le ighton i quickly suggested that the 

random-walk of magnetic elements by the evolv- 
ing supergranules could transport following polar- 
ity elements poleward to reverse the Sun's polar 
field each solar cycle while leading polarity ele- 
ments would be transported across the equator 
where they would meet and cancel their opposite 
polarity counterparts. While we now recognize 
that Leighton over estimated the effective diffu- 
sivity implied by this random walk (and that the 
Sun's poleward meridional flow provides the miss- 
ing transport), supergranule diffusion is still not 
fully constrained or understood. Leighton's initial 
estimate of 770 — 1540 km 2 s -1 for the diffu sivity 
was reduced to 200 — 400 km 2 s _1 by iMosher 



(|19771 ) who studied the dis placements of chro 



mosph eric network elements. ISchriiver fe Martin 
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(199fJ) measured the displacements of the mag- 
netic elements themselves using magnetograms 
and found a value of 250 km 2 s _1 in plage sur- 
rounding active regions but only 110 km 2 s _1 in 
the active regions themselves. Models for the 
surface magnetic flux tra nsport typically use val- 
ues ab out twice this size - Ivan Ballegooiien et al 
(|l998l ) used 450 km 2 s" 1 while lWang et al. I (| 200 
used 500 km 2 s _1 . 

Supergranules are embedded in a surface shear 
l ayer of t heir own maki ng. Foukal fe Jokipiil 
(|l975h and iFoukal I (1 19791) suggested that the ra- 
dial flows within the supergranules would con- 
serve angular momentum and produce a shear 
layer with slower rotation at the surface and more 
rapid rotation bel ow - consistent wi t h the ear- 
lier observations of Howard fe Harvey ( 1970t) that 
the sunspots rotate about 5% faster than the sur- 
rounding surface plasma. The production of this 



shear lay er was succesfully modell ed in spherical 
shells by iGilman fc Foukal 1 (Il979l) and in plane- 
parallel layers bv lHathawav I (jl982h . The theo- 
retical arguements suggest a change in rotation 
frequency such th at 8fl ~ —2Q,6r/r. In addition, 
lHathawav I (jl982h found that a meridional flow 
was also produced with poleward flow near the 
surface and equatorward flow at the bottom. 

The surface shear layer has been probed with 
helioscismic techniques which confirm the exis- 
tence of the layer but give somewhat inconsistent 
results for the detailed structure in depth and lat- 
itude. Thompson et al. ( 19961 ) used the global 
oscillation modes from GONG and found that the 
rotation rate increased inward to a depth of about 
35 Mm (Sr/r ~ 5%) and that the rate of in- 
crease was larger at t he eq uator than at 30° lat- 
itude. ISchou et al. I |l998) found similar results 
from global modes in MDI data but noted that 
the shear appeared to disappear by 60° latitude 
and perhaps reverse sign at higher latitudes. 

Basu et al. (1999) used the local helioseismol- 
ogy technique of ring-diagram analysis. They 
found a thin (4 Mm) outer shear layer in which the 
shear did not appear to reverse at high latitudes 
and a deeper shear layer in which it did. They also 
found a poleward meridional flow which increased 
in amplitude across the thin outer shear layer but 
then remained nearly constant across the inner 
shear layer. Their meridional flow had peak veloci 



ties of about 50 m s at 50° latitude. Giles et al 



(|l997l) used a time-distance local helioseismolog 



y 

technique and found a meridional flow with a peak 
velocity of only about 24 ms" 1 at 45° latitude and 
indicated that it was constant to a depth of about 
26 Mm. 

The observed rotation of the supergranules has 
added further mystery to their nature. (The mech- 
anism that selects the characteristic size o f super 
granu les is perhaps the biggest mystery.) iDuvall 



(1980) cross-correlated the equatorial Doppler ve- 
locity patterns and found that the supergran- 
ules rotated more rapidly than the plasma at the 
photosphere and that even faster rotation rates 
were obtained when there were longer (24-hour 
vs. 8-hour) time intervals between Doppler im- 
ages. He attributed this behavior to the surface 
shear layer in which larger, longer-lived, cells ex- 
tend de^£e^Jnto_rtie_jnore rapidly rotating lay- 
ers. Snodgrass fc Ulrich I (|l990l ) used data from 
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Mount Wilson Observatory to find the rotation 
rate at different latitudes and noted that the rota- 
tion rates for the Doppler pattern were some 4% 
faster than the spectroscopic rate and, mysteri- 
ously, some 2% faster than the small magnetic 
features that are o bserved to outl i ne th eir bor- 
ders. More rece ntlv lBeck fc Schou I (l200flh used a 
Fourier transform method to find that the larger 
features do rotate more rapidly than the smaller 
features and that the low wavenumber components 
of the Doppler pattern rotate more rapidly than 
the plasma at any level within the surface shear 
layer. This led them to suggest that supergranules 
have wave-like characteristics with a preference for 
prograde propagation. 

Hathaway et al. ( 20061 ) showed that this "super- 



rotation" of the Doppler pattern could be at- 
tributed to projection effects associated with the 
Doppler signal itself. As the velocity pattern ro- 
tates across the field of view the observed line-of- 
sight component is modulated in a way that essen- 
tially adds another half wave to the pattern and 
gives a higher rotation rate that increases propor- 
tionally with decreasing wavenumber. They took 
a fixed velocity pattern (which had spatial charac- 
teristics that matched the SOHO/MDI data) and 
rotated it rigidly to show this "super-rotation" 
effect. While this indicated that the Doppler pro- 
jection effect should be accounted for, the fixed 
pattern coul d not account for al l the variations 
reported bv lBeck fe Schou I (|2000h . 

Gizon et al. I ([20031 ) used time-distance helio- 



seismology to find the supergranular flow (rather 



than direct Doppler measurements) and ISchou 
((2003J) "divided-out" the line-of-sight modula- 
tion. Both studies found slower rotation rates 
that matched that of the magnetic features but 
saw evidence for wave-like prograde and retro- 
grade moving components. iMeunier fc Roudier 



(|2007l ) compared rotation rates obtained by track- 
ing Doppler features, magnetic features, and diver- 
gence features (which were, in turn, derived from 
correlation tracking of smaller intensity features) . 
While they concluded that projection effects do 
influence the rotation rate determined from the 
Doppler features, they found that the magnetic 
features rotate more slowly than the supergran- 
ules. 

In this paper we report on our analyses of data 
from the SOHO/MDI instrument Scherrer et al 



(|l995l )] and from simulated data in which the su- 
pergranules are simply advected by differential ro- 
tation and meridional flow that vary with latitude 
and depth. The simulated data are designed to 
faithfully mimic the SOH O/MDI data [th e same 
dat a that was an alyzed in Beck fc Schou ( 2000h 
and lSchou I (|2003h ] with simple assumptions about 
the dynamical structure of the surface shear layer. 
The analyses include reproductions of those done 
in earlier studies. Through the simulations we can 
better determine the actual differential rotation 
and meridional flow profiles consistent with the 
Doppler observations (which are subject to line- 
of-sight projection effects). 

2. The Data 

The full-disk Doppler images from SOHO/MDI 
are obtained at a I-minute cadence to resolve the 
temporal variations ass ociated with the p -mode 
oscil lations. We [c f. Hathaway et al. ( 2000l) 



and Beck fc Schou ( 200ol )] have temporally fil- 



tered the images to remove the p-mode signal by 
using a 3I-minute long tapered Gaussian with a 
FWHM of 16 minutes on sets of 31 images that 
were de-rotated to register each to the central im- 
age. These filtered images were formed at 15- 
minute intervals over the 60-day MDI Dynamics 
Run in 1996. This filtering process effectively re- 
moves the p-mode signal and leaves behind the 
Doppler signal from flows with temporal variations 
longer than about 16 minutes. 

Supergranules, with typical wavenumbers of 
about 110, are very well resolved in this data (at 
disk center wavenumbers up to 1500 are resolved). 
While granules are not well resolved, they do ap- 
pear in the data as pixel-to-pixel and image-to- 
image "noise," as a convective blue shift (due to 
the correlation between brightness and updrafts), 
and as resolved structures for the largest mem- 
bers. These data are prepared for studying the 
cellular features by first measuring and remov- 
ing the Doppler signals due to: 1) observer mo- 
tion, 2) convective blue shift, 3) differential rota- 
tion, and 4) the axisymmetricmeridional flow. The 
data are then mapped onto heliographic coordi- 
nates with equal spacing in both longitude and lat- 
itude. This mapping includes accounting for the 
position angle and tilt angle of the Sun's rotation 
axis as well as the orientation of the SOHO/MDI 
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detector. Note that the position angle and tilt 
angle of the Sun' s rotation axis determined by 
Carrington I (|l863l ) have been found to b e in er- 
ror by about 0.1° bvlBeck fc Giles I ([20051) and by 



Hathaway fc Rightmire ( 2010h . This correction is 



included in this study. 

The simulated data are constructed in the man- 
ner described previous l y (lHathaway" 19881 1992 ; 
Hathaway et al. 200Ct l2002h from vector veloci- 
ties with 
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where YJ n (0, cj>) is the spherical harmonic function 
of degree I and azimuthal order m , 9 is the colati- 
tude measured southward from the north pole, and 
<f> is the longitude measured prograde from the cen- 
tral meridian or some fixed longitude. The com- 
plex quantities R™, 5™, and T™ are the spectral 
coefficients for the radial, poloidal, and toroidal 
components, respectively. 

To simulate the observed line-of-sight velocity, 
the spectral coefficients from I = to £ — 1500 are 
specified and the three vector velocity components 
are calculated using Eqns. 1-3 on a grid with 1500 
points in latitude and 4096 points in longitude. 
The Doppler velocity image is constructed by de- 
termining the longitude and latitude at a point on 
the image, finding the vector velocity at that point 
using bi-cubic interpolation, and then projecting 
that vector velocity onto the line-of-sight using 



V, 



los \ 



V r (9, <f>) sin B cos( 



V r (9,(j)) cos B sin# cos^ 
Ve (0,4>) sin B sin 6 - 
Ve (9, 4>) cos Bq cos 9 cos i 
Vj,(6, 4>) cos i?o sin0 



(4) 



where Bq is the latitude at disk center (or equiv- 
alently the tilt of the Sun's north pole toward the 
observer) and velocities away from the observer 
are taken to be positive. The line-of-sight velocity 
at an array of 49 points within each pixel is de- 
termined and an average is taken to simulate the 
integration over a pixel with the MDI instrument. 

With the current simulations the instrumental 
blurring is treated in a more realistic manner. Pre- 
viously the Doppler velocity image itself was con- 
volved with an MDI point-sprcad-function. Here 
we make red and blue intensity images from our 
Doppler velocity image and a simple limb dark- 
ened intensity image, convolve those with an MDI 
point-spread-function, and construct the blurred 
Doppler velocity image from the difference divided 
by the sum. 

This process yields Doppler velocity images 
that are visually indistinguishable from MDI 
Doppler velocity images. A velocity image com- 
parison is shown in Fig. 1. 



Fig. 1. — Doppler velocity images (MDI on left, 
simulation on right) in which the observer's mo- 
tion, the convective blue shift, and the Doppler 
signals due to the axisymmetric flows — differen- 
tial rotation and meridional circulation — are re- 
moved to better reveal the cellular structures (su- 
pergranules). The spatial characteristics are visu- 
ally indistinguishable. 

The velocity pattern is evolved in time by intro- 
ducing changes to the spectral coefficients based 
on two processes - the advection by the axisym- 
metric flows (differential rotation and meridional 
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flow) and random processes that lead to the finite 
lifetimes of the cells. 

The advection is governed by an advection 
equation 



dw 
~dt 



V{9)dw U{9) dw 



89 



r sin ( 



(5) 



where w is a velocity component, U(9) = r sin 9Vt{9) 
gives the differential rotation profile and V(9) 
gives the meridional flow velocity profile. Repre- 
senting w as a series of spherical harmonic com- 
ponents (Eqs. 1-3) and projecting this advection 
equation onto a single spherical harmonic gives 
a series of coupled equations for the evolution of 
the spectral coefficients (Appendix A) . Solid body 
rotation simply introduces a constantly increasing 
phase for each coefficient. Differential rotation 
couples the phase change in one spectral coeffi- 
cient to spectral coefficients with wavenumbers 
£ ± 2 and £ ± 4 for differential rotation of the form 



n{9) = n + n 2 cos 2 9 + n 4 cos 4 e 



(6) 



while a simple but reasonable meridional flow pro- 
file with 



V{9) = V cos9sm9 



(7) 



couples one spectral coefficient to spectral coeffi- 
cients with wavenumbers £±2. (Spherical harmon- 
ics have fixed latitudinal structure. Spectral power 
must pass from one spherical harmonic component 
to another in order to move a feature in latitude.) 

These cellular flows are embedded in the Sun's 
surface shear layer. We approximate the change in 
the rotation ra te in the out e rmost 5% of the Sun 
as reported bv lHowe et al. ( 2007 ) with 



fi(r,0) = fl(0)/(r) 



where 



f(r) = 1 + 0.038 



1 - e 
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and the latitude dependence is given by 
— = 452 nHz (14.07° day" 1 ) 
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(11) 



(12) 



Assuming that the cells extend to depths similar 
to their horizontal dimensions, and that they are 
advected at flow rates representative of that depth, 
Eq. 9 is transformed into a function of £ with 



f{£) = 1 + 0.038 



1 _ e -°*H/e 



(13) 



This shear layer profile is illustrated in Fig. 2 
along with the gradients expected from theoret- 
ical arguements for flows that conserve angular 
momentum. We assume a meridional flow which 
is constant with depth across this layer and has a 
latitude dependence characterized by 



-30 m 



(14) 



which gives a peak meridional flow velocity of 
15 m s" 1 at 45° latitude. 

The finite lifetimes for the cells are simulated 
by introducing random perturbations to the spec- 
tral coefficient phases. The size of these perturba- 
tions increases with wavenumber to give shorter 
lifetimes to smaller cells with 



6$T oc \/At/r(£) 



(15) 



where 5$™ is the change in phase for a complex 
spectral coefficient of degree £ and order m, At is 
the time interval between simulated Doppler im- 
ages, and t(£) is proportional to the lifetime for 
a spectral component of degree £. Lifetimes are 
well approximated by a turn-over time for turbu- 
lent convective flows. The cellular flow velocities 
are roughly proportional to £ while their diame- 
ters are inversely proportional to £. The turn-over 
times should then be inversely proportional to £ 2 . 
We find a reasonable fit to the data using 



T [l) = 6.5^hrs 



P 



(16) 



3. The Analyses 



Several anaylsis programs were applied to both 
the MDI data and the simulated data. Power spec- 
tra were obtained to characterize and compare the 
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Fig. 2. — The assumed rotation rate as a func- 
tion of radius in the surface shear layer for three 
different latitudes. The surface rotation rate at 
each latitude is indicated by the dashed lines. The 
variations in rotation rate for flows that conserve 
angular momentum from the surface inward are 
indicated by the dotted lines. 

distribution of cell sizes and flow velocities. The 
rotation of the Doppler pattern was determined 
using multiple techniques based on previous stud- 
ies. The meridional flow was measured based on 
the movement of the Doppler pattern and cell life- 
times were estimated from cross-correlation anal- 
yses. 

3.1. Convection Spectra 

Convection spectra for individual im ages were 
obtain ed us ing the methods de scribed bv lHathawav I 
(|l987l ) and lHathawav I (| 19921) - the Doppler sig- 
nal due to the motion of the observer is removed, 
the convective blue shift signal is identified and 
removed, the data is mapped to heliographic co- 
ordinates, the axisymmetric flow signals due to 



differential rotation and meridional circulation are 
identified and removed, and the remaining signal 
is projected onto spherical harmonics. The aver- 
aged spectra from the 60-day 1996 MDI Dynamics 
Run and from our 10-day simulated data run are 
shown in Fig. 3. 

The match between these spectra is obtained 
by adjusting the input spectrum for the simulated 
data. This spectrum contains two Lorentzian-likc 
spectral components - a supergranule component 
centered on £ ~ 110 with a width of about 100 
and a granule component centered on £ ~ 4000 
with a width of about 3000. The MDI spectrum 
is well matched with just these two components 
without the additi on of the mesogra n ule co mpo- 



nent suggested by iNovember et al. ' (1981). In 



fact, we have a distinct dip in the input spectrum 
at wavenumbers £ ~ 500 that should be re p resen - 
tative of mesogranules. Hathaway et al. ( 2000h 



showed that this dip is easily seen in the observed 
spectra from the MDI high resolution data. 
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Fig. 3. — Power spectra from spherical harmonic 
analyses of MDI data (solid line) and simulated 
data (dotted line). The spectral match indicates 
that both datasets have very similar spatial struc- 
tures. (The rapid drop in power for wavenumbers 
beyond 600 is largely due to instrumental blur- 
ring.) 

3.2. Rotation 

Rotation information is obtained when addi- 
tional analyses are applied to the data after it 
has been mapped onto heliographic coordinates. 
Longitudinal strips of this data, centered on lat- 
itudes from 75° south to 75° north, were cross- 
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correlated with correspo nding; st r ips fr om later 
images as was done by [ Duvall (1980) and by 



Snodgrass fc Ulrich I (|l99dh . The shift in the lo- 
cation of the cross-correlation peak divided by the 
time difference gives the rotation rate. 

Fig. 4 shows the differential rotation from the 
cross-correlation analysis. The profiles from the 
simulated data match those from the MDI data at 
all but the highest latitudes. Both datasets show 
fast er rotat io n rate s for long er time lags as noted 
bv iDuvall I (|l980f ) and by ISnodgrass fc Ulrich 
(1990). The match between the MDI data and 



the simulated data indicates that the latitudi- 
nal differential rotation profile is fairly well rep- 
resented by Eq. 6 with the coefficients given by 
E qns. 10-12 derived from the helioseismic studies 
of Howe et al. ( 20071 ) . (The simulated data would 
be a better match at high latitudes for H4 — —85 
nHz.) The MDI data does show slightly faster 
equatorial rotation. This indicates that either Eq. 
13 should be modified to give faster rotation at the 
lower wavenumbers or Eq. 16 should be modified 
to give shorter lifetimes at the higher wavenum- 
bers. 
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Fig. 4. — Differential rotation profiles from cross- 
correlation analyses of MDI data (solid lines) and 
simulated data (dotted lines) for four different 
time lags between Doppler images - 2, 4, 8, and 
16 hours. All profiles are well matched at all time 
lags and at all but the highest latitudes. 



We reproduced the analysis of iBeck fc Schou 
(2000). Lines of data from the heliographic maps 



were Fourier analyzed in longitude and those spec- 
tral coefficents were then Fourier anaylzed in time 
over 10-day intervals. The power spectra were av- 



eraged over all latitude lines between ±9°. The 
rotation rate for each wavenumber t was deter- 
mined by finding the temporal frequency of the 
center of gravity of the power using a frequency 
window of ~ 58/xHz (and iterating on the position 
of that window four times) and then dividing the 
temporal frequency by the wavenumber. 

Fig. 5 shows the equatorial rotation rates as 
functions of wavenumber for both the MDI data 
and the simulated data. Here again we find a 
good match for all but the lowest wavenumbers. 
This indicates that Eq. 13 gives a good wavenum- 
ber dependence for the rotation rate. The drop 
in rotation rate for the MDI data at wavenum- 
bers below 30 is due to imaging artifacts that in- 
troduce signal at low temporal frequencies which 
can overlap with the rotation signal at low spatial 
frequencies. The effects of the line-of-sight pro- 
jec tion on measured rotat ion rates, as discussed 
bv lHathawav et al. I ([20061) . are evident in the in- 
crease in the observed rotation rates above the 
input rates from Eq. 13 for wavenumbers below 
~ 100. Although the higher noise level from the 
single 10-day simulation makes precise compar- 
isons difficult, it does appear that the MDI data 
might be better fit if the rotation rates given by 
Eq. 13 were somewhat higher at wavenumbers be- 
low - 100. 
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Fig. 5. — Equatorial rotation rates as functions 
of wavenumber from analyses of 6 10-day sets of 
MDI data (solid line) and a single 10-day set of 
simulated data (dotted line) . The equatorial rota- 
tion rates used as input (Eq. 13) are shown with 
the dashed line. The rotation rates match very 
well at all but the lowest wavenumbers. Rotation 
rates faster than the input rates at wavenumbers 
below about 100 are due to projection effects. 
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We hav e also reproduced one of the analyses 
of ISchou I (|2003r ) . Lines of data from the helio- 
graphic maps are apodized and then multiplied by 
longitude dependent weighting functions designed 
to remove the Doppler projection effect and to 
isolate either longitudinal motions or latitudinal 
motions. The weighting used to isolate longitu- 
dinal flows is based on the final term in Eq. 4 
but with B = and a constant of 0.01 added to 
sin <fi to avoid division by zero at the central merid- 
ian. The lines are shifted in longitude according 
to the differential rotation rate and then Fourier 
analyse d in space and time to obtain "fc — u" di- 



agrams. ISchou ( 2003 ) noted that these diagrams 
show power for both prograde (ui/k < 0) and ret- 
rograde (oj/k > 0) motion but with excess power 
in the prograde components. 

Fig. 6 shows the k — u> diagrams from 10 days 
of the MDI data and from the simulated data for 
latitudes between ±4.5°. Both show prograde and 
retrograde components as similarly shaped wedges 
of enhanced power extending from the origin to 
both positive and negative longitudinal wavenum- 
bers. However, the MDI data does show excess 
power for prograde motion (power excess extend- 
ing from the upper left to the lower right) that is 
not evident in the simulated data. (Note that run- 
ning this analysis on the full 60-day run of MDI 
data shows a clear separation between prograde 
and retrograde motions.) 

3.3. Meridional Flow 

The meridional motion of the Doppler pat- 
tern can be measured by the cross-correlation 
method used to find the differential rotation pro- 
files illustrated in Fig. 4 or by usi ng the apodiz - 
ing/ weighting method described bv lSchoul ( 2003h . 
Fig. 7 shows a comparison of the meridional flow 
profiles obtained from the MDI data and the sim- 
ulated data when we take Vq = —30 m s _1 , repre- 
senting a poleward meridional flow with peak ve- 
locity of 15 m s _1 . The meridional flow we find for 
the Doppler features in the simulation very closely 
matches the input meridional flow (dashed lines in 
Fig. 7). This indicates that the projection effects 
that produce the super-rotation of the Doppler 
features do not impact the measurements of the 
meridional flow. 

It is clear that the actual meridional flow of 
the supergranules is more complicated than that 
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Fig. 6. — The k — u> diagrams from analyses 
of MDI data (left) and simulated data (right). 
Darker shades indicate more power. Both datasets 
show prograde and retrograde components with 
similar power and limits but the MDI data shows 
excess power in the prograde components (power 
in upper left and lower right quadrants) . The thin 
line from upper right to lower left in the MDI data 
is produced by imaging artifacts which appear to 
move retrograde at the Sun's rotation rate. 
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Fig. 7. — Meridional flow profiles from cross- 
correlation analyses of MDI data (solid lines) and 
simulated data (dotted lines) for four different 
time lags between Doppler images - 2, 4, 8, and 
16 hours. The meridional flow profile input to the 
simulation is shown by the dashed lines. 



used in our simulation. The meridional flow pro- 
files for the supergranules shown in Fig. 7 for the 
MDI data ag r ee ver y well wi th those obtai ned by 



Gizon et al. 



(|2003l ) and by ISchou I (|2003l ) using 
somewhat different methods. The MDI data in- 
dicates the presence of counter-cells in the polar 
regions during the 1996 data run that would be 
better represented by 
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V(ff) — (50 cos 9 — 60 cos^ 9) sin 9 (17) 

In addition the flow appears to weaken for longer 
time differences between the cross-correlated data. 
This suggests the the meridional flow in the sur- 
face shear layer decreases in amplitude with depth. 

3.4. Lifetimes 

The lifetimes of the cells can be estimated by 
comparing the strength of the cross-correlation co- 
efficients for the Doppler features as functions of 
the time interval between Doppler images. Fig. 8 
shows the strength of the correlations as functions 
of latitude and time lag for both the MDI data 
and the simulated data. The profiles are similar 
in shape but slightly flatter with the MDI data. 
The correlation coefficients are well matched at 
the equator for time differences of 4-hours and 
16- hours but the simulation shows stronger cor- 
relations at 2-hours and weaker correlations at 
8-hours. This suggests that the formulae (Eqns. 
15 and 16) for the random phase changes to the 
spectral coefficients needs slight modification. In 
particular the cellular features that dominate at 
At = 2 hours need shorter lifetimes. 
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Fig. 8. — Cross-correlation peak profiles from 
analyses of MDI data (solid lines) and simulated 
data (dotted lines) . The profiles are well matched 
at At — 4 and At — 16 hours but with slight 
differences at At — 2 and At = 8 hours. 



4. Conclusions 

We have produced simulated data in which the 
cellular structures (supergranules) are advected 



by differential rotation and meridional flow and 
evolved by uncorrelated random changes. When 
we compare results from analyses of these data 
with those from the MDI data we find that the 
simulated data exhibit much of the same char- 
acteristics as the MDI data - the visual struc- 
tures, the power spectra, the rotation, the merid- 
ional flow, and the evolution rates all match fairly 
closely. While some of the rotation characteris- 
tics have_j2eeii_j^ributed t o wave-like pr o pertie s 
by iBeck fc Schoul (<2000h . iGizon et al. 1 d2003h . 



and ISchou I ( |2003f ), our simulated data is simply 
advected by a zonal flow (differential rotation) 
with speeds that never exceed thos e determined 
from helioseismology as reported by Irlowe et al. I 
<|2007h . 

While the similarities between the simulated 
data and the MDI data are strong, the differ- 
ences are important and interesting. The rota- 
tion variations with depth and latitude were cho- 
sen to closely mat ch those shown in Fig 1. of 
Howe et al. I ( 2007 ). The differences in rotation 
rate seen at high latitudes in our Fig. 4 sug- 
gest that the magnitude of SI 4 should be increased 
from -75 nHz to -85 nHz. The differences in rota- 
tion rate seen at low wavenumbers in our Fig. 5 
suggest an increase in rotation across the surface 
shear layer of about 2 nHz at the equator. 

The lack of excess power in the prograde com- 
ponents represented in Fig 6 for the simulation 
and the lack of evidence for a split between the 
prograde and retrograde components indicate dif- 
ferences in the evolution of the pattern. The man- 
ner in which we evolve the cells in the simulation 
is not fully satisfactory. The phase (and ampli- 
tude) changes associated with the evolving cells 
should replicate their advection by large r (non - 



axysymmetric) flows ( Williams fc Cuntz 2009), 



their break-up into smaller cells, and their join- 
ing to form larger cells. The random changes in 
phase that we introduce to evolve the cells pro- 
duces random offsets that are larger for larger 
wavenumbers and, when modulated by the power 
spectrum shown in Fig 3, produces the wedge of 
power shown in Fig 6 for the simulation. The split 
between prograde and retrograde movement seen 
in the MDI data suggests that the mo tions are not 
rando m but preferentially east- west. iLisle et al. 



(|2004h found that supergranules appear to have 
a persistent north-south alignment and suggested 
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that this was due to a north-south elongation of 
larger, giant cells that advect the supergranules 
to their boundaries. The rotationa l const raints on 
these giant cells [cf. iMiesch et al. (2008) and ref- 
erences therein] should in-fact give a preference 
for east-west flows near the equator and may pro- 
duce an added preference for prograde flows. Note 
that this interpretation of the results does not in- 
volve wave-like properties for supergranules - but 
instead invokes simple advection by larger flows. 

The advection of the supergranules by the 
meridional flow shows considerable promise. The 
fact that we fully and accurately recover the input 
velocity profile indicates that the latitudinal move- 
ment of supergranules provides new information 
on the meridional flow. The results with the MDI 
data indicate the presence of counter-cells (equa- 
torward flow) at latitudes above about 65° — 70°. 
These results also indicate a meridional flow speed 
that matches that found for the small magnet ic el- 
em ents as measured by Komm et al~l (jl993l ) and 
by lHathawav k. Rightmire ( 2010 ) but is signif- 
icantly slower than that reported in the helio- 
seismic dete r minat ion of iGiles et al. I (|1997l ) and 
Basu et al. The indication of a decrease 

in meridional flow speed with depth is also at odds 
with those helioseismic results but in agreement 



with the simulations of Hathaway (198 



Magnetic element positions are known to co- 
incide with the boundaries of supergranules and 
move as the supergranule bou ndaries evolve. Thi s 
was shown fairly explicitly by iLisle et al~ ] (<2000h . 



This suggests that the advection of the super- 
granules by the axisymmetric flows should be di- 
rectly related to the advection of the magnetic el- 
ements. The differential rotation itself (relative to 
the Carrington rotation) has a velocity range from 
~ 40 m s _1 prograde to ~ 180 m s _1 retrograde 
while the meridional flow has a peak velocity of 
only ~ 11 m s -1 . These weak flow velocities are 
no match for the flows of several 100 m s _1 in the 
supergranules. The magnetic elements should be 
quickly carried to the boundaries and ultimately 
the interstices of the supergranules where the di- 
rect action of the differential rotation and merid- 
ional flow are too weak to dislodge them. The 
magnetic elements should experience differential 
rotation and meridional flow only to the extent 
that the supergranules themselves are advected by 
these flows at velocities representitive of deeper 



layers. 

The variation in the differential rotation and 
meridional flow with depth and latitude within 
the surface shear layer needs further examina- 
tion. We have already argued for some changes 
from what was assumed in our simulation. In ad- 
dition, we note that our association of different 
depths with different wavenumbers £ as indicated 
in Eqns. 9 and 13 suggests that cells with diame- 
ters D = 2ttRq/£ are advected by flows at depth 
i?© — r 0.82D. This would make typical super- 
granules extend to depths of 30 Mm. While this is 
deeper than suggested by t he he lioseis mic inves- 



tigatio ns of iDuvall et al. I (|1997j ) and IWoodard I 
(|20070 , those investigations were limited to lay- 
ers not much deeper than 7 Mm. Shallower cells 
would suggest that the surface shear layer is shal- 
lower as well. 
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A. Spectral Coefficient Changes Due to Advection 

The changes in the spectral coefficients produced in the advection of the pattern by the axisymmetric 
flows are best illustrated using the radial component of the flow pattern as given by 



■;(M,t)= E E R?(t)Pp(x)e 



im<fi 



(Al) 



=0 m=-t 



where the time- varying spectral coefficient is R™(t), x = cos 6, and P™ is an Associated Legendre polynomial 
of angular degree I and azimuthal order to. The Associated Legendre polynomials are normalized such that 



(A2) 



PT(x)PF (x)dx = 6} S% 

The advection of this flow pattern by the axisymmetric flows is given by 

dw V dw U dw 
dt r d9 rs'mO d<fi 

The spherical harmonic representation gives explicit expressions for the spatial derivatives with 



q Imax £ ^ 

i§ = E E R w^ { 



=0 m=-t 



H+l 



im(p 



and 



dw 
~d4) 



Imax l 



where 



Taking the meridional velocity as 



and the azimuthal velocity as 



E E RTWrnPrixy 

e=o m=-l 

m= r (2^ + i)(2f - 1) 

[{£ + m)(£-m) 
V{6) = V Q cos 9 sin 6 



im(j> 



U(9) = r sin 6>(O + Vl 2 cos 2 6 + n 4 cos 4 9) 
and then projecting eq. A3 onto YJP 1 gives 



dt 



T/ r Imax p\ 

e=m' J ~ 1 

Imax „i 

-im'n R™' - im'fl 2 E R ™' / x 2 pf {x)pf (x)d2 

Imax -.1 

-im'n 4 J2 R T'j x 4 Pf {x)Pf {x)dx 



dx 



(A3) 



(A4) 



(A5) 



(A6) 



(A7) 



(A8) 



(A9) 
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The products of x, x 2 , and x 4 with the Associated Legendre polynomials can be replaced with expressions 
containing only Associated Legendre polynomials of the same azimuthal order, m, using the recursion relation 



PP(x) = AfxP^x) - BfP?_ 2 {x) 



(A10) 



where A™ is given by eq. A6 and 



B't 



{2£+\){£ + m- l)(£-m- 1) 



{2£ - 3)(£ + m)(£ - m) 
Eq. A9 then reduces to (after dropping the primes) 



(All) 



9RJ 
dt 



Vo 
r 



(£-2) 



ATA 



RYL-2 + 



-ZTOf2 4 

-imCl4 



1 



IB 



£+2 



Am Am 
A £+2 A £+1 



ATA 



-R 



t-2 



+ 



1 



B 



1+2 



Am Am 
A £ A £ 

BV 



It f 



Am Am 

Ji e+2 Ji e+2 



rt, 



t+2 



Am Am 
A £+2 A £+1 



+ 



Am An 
A l+l A l 



RT + 



B 



e+3 
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Ji £+3 Ji £+2 
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D £+2 IX £+2 
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Am An 
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B 



on 

■SXn 
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+ 
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(A12) 



A common factor in Eq. A12 is of the form 

BT 



!)(£■ 



1) 



ATA 



e-i 



(2£- l)(2£-3) 



The dominant terms on the RHS of Eq. A12 are those for the direct rotation of the RT coefficients, 
terms are best handled analytically by taking 



(A13) 
These 



RT = TiTe im ^ 1 



(A14) 



where 



9, r i 



fl 4 



^2 

B 



B. 



£+2 



B. 



£+1 



Am Am 
A £+2 A l+\ 
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A £+l A £ 
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B: 
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+ 
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ATA 



£ ^£-1 



B 
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(A15) 
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With this substitution, Eq. A12 becomes an equation for the evolution of 1Z™ in which the meridional 
flow (terms multiplied by Vb) couples the amplitude of the £ component to the £ ± 2 components while the 
differential rotation (terms multiplied by O2 and CI4,) couples the phase of the £ component to the £ ± 2 and 
I ± 4 components. These spectral coefficients can be evolved in time using Eq. A12 with a fourth-order 
Runge-Kutta scheme. The solutions are stable and well behaved for quarter-hour time steps over at least 
10-days for £ upto at least 1500. 
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